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Abstract The Henon-Heiles potential is undoubtedly one 
of the most simple, classical and characteristic Hamiltonian 
systems. The aim of this work is to reveal the influence of 
the value of the total orbital energy, which is the only pa¬ 
rameter of the system, on the different families of orbits, by 
monitoring how the percentage of chaotic orbits, as well as 
the percentages of orbits composing the main regular fam¬ 
ilies evolve when energy varies. In particular, we conduct 
a thorough numerical investigation distinguishing between 
ordered and chaotic orbits, considering only bounded mo¬ 
tion for several energy levels. The smaller alignment index 
(SALI) was computed by numerically integrating the equa¬ 
tions of motion as well as the variational equations to exten¬ 
sive samples of orbits in order to distinguish safely between 
ordered and chaotic motion. In addition, a method based on 
the concept of spectral dynamics that utilizes the Fourier 
transform of the time series of each coordinate is used to 
identify the various families of regular orbits and also to rec¬ 
ognize the secondary resonances that bifurcate from them. 
Our exploration takes place both in the physical (x, y) and 
the phase (y, y) space for a better understanding of the orbital 
properties of the system. It was found, that for low energy 
levels the motion is entirely regular being the box orbits the 
most populated family, while as the value of the energy in¬ 
creases chaos and several resonant families appear. We also 
observed, that the vast majority of the resonant orbits belong 
in fact in bifurcated families of the main 1:1 resonant fam¬ 
ily. We have also compared our results with previous similar 
outcomes obtained using different chaos indicators. 
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1 Introduction 

The interest in the existence of a third integral of motion 
for stars moving in the potential of a galaxy had been re¬ 
vived back in late 50s and in the early 60s. It was initially 
assumed that the potential had a symmetry and it was time- 
independent, therefore in cylindrical coordinates (r, 6, z) it 
would be only a function of r and z. There should exist 
five integrals of motion that are constant for six-dimensional 
phase space. However, the integrals can be either isolating or 
non-isolating. The non-isolating integrals normally fill the 
entire available phase space and give no restriction to the 
orbit. 

By the time when Henon and Heiles wrote their pio¬ 
neer paper, there were only two known integrals of motion 
namely the total orbital energy and the angular momentum 
per unit mass of the star. It can easily be shown, that at least 
two of the integrals, in general, are non-isolating. It was also 
assumed, that the third integral was also non-isolating be¬ 
cause no analytical solution had been found for it so far. 
Nevertheless, observations of stars near the Sun as well as 
numerical computations of orbits in some cases behaved as 
if they obeyed in three isolating integrals of motion. 

Henon and Heiles set out to see if they could find any 
real proof of that there should exist a third isolating integral 
of motion. They did this by conducting numerical computa¬ 
tions, even though they did not hold too hard to the astro¬ 
nomical meaning of the problem; they only demanded that 
the potential they investigated was axial symmetric. They 
also assumed that the motion was confined to a plane and 
went over to Cartesian phase space (x,y,x,y). After some 
trials, they managed to come up with a valid potential. This 
potential is analytically simple so that the orbits could be 
calculated rather easily but it is still complicated enough so 
that the types of orbits are non trivial. This potential is now 
known as the Henon and Heiles potential ED. 
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Over the years, a huge load of research has been devoted 
on the Henon-Heiles system (e.g., |[38l[T7l[23l[29l l9imi^[3l 
|6l[ia|30l[3ll[7l|3211[8l|33[3l|^ At this point we 
should emphasize, that all the above-mentioned references 
regarding previous studies in the Henon-Heiles system are 
exemplary rather than exhaustive, taking into account that a 
vast quantity of related literature exists. The vast majority 
of these papers deals mostly either with the discrimination 
between regular and chaotic motion or with the escape prop¬ 
erties of orbits. In the present paper, on the other hand, we 
proceed one step further contributing to this active field by 
classifying ordered orbits into different regular families and 
also monitoring how the percentages of the main families of 
orbits evolve as a function of the total orbital energy, con¬ 
sidering only energy levels below the escape energy, where 
orbits are bounded and a variety of chaotic and several types 
of ordered motions exist. This is exactly what makes this 
paper novel. 

In a thorough pioneer study 1^ . analyzed the orbital 
content in the coordinate planes of triaxial potentials but also 
in the meridional plane of axially symmetric model poten¬ 
tials, focusing on the regular families. Few years later, ca 
developed a method based on the analysis of the Fourier 
spectrum of the orbits which can distinguish not only be¬ 
tween regular and chaotic orbits, but also between loop, box, 
and other resonant orbits either in two or three dimensional 
potentials. This spectral method was improved and applied 
in Eli in order to identify the different kinds of regular or¬ 
bits in a self-consistent triaxial model. The same code was 
improved even further in B2l . when the influence of the cen¬ 
tral nucleus and of the isolated integrals of motion (angular 
momentum and energy) on the percentages of orbits in the 
meridional plane of an axisymmetric galactic model com¬ 
posed of a disk and a spherical nucleus were investigated. 

The structure of the present paper is as follows: In Sec¬ 
tion we present a detailed description of the properties of 
the Henon-Heiles system. All the different computational 
methods used in order to determine the character as well 
as the classification of the orbits are described in Section 
In the following Section, we conduct a thorough analy¬ 
sis of several sets of initial conditions of orbits presenting 
in detail all the numerical results of our computations. Our 
article ends with Section where the discussion and the 
conclusions of this research are presented. 


2 Properties of the Henon-Heiles system 

One of the most paradigmatic model potentials for time- 
independent Hamiltonian systems with two degrees of free¬ 
dom (2 d.o.f) is undoubtedly the Henon-Heiles system ll24l . 
This particular system defines the motion of a test particle 


with unit mass and the corresponding two-dimensional po¬ 
tential is given by 

V(x,y)=^{x^+y^) + x^y-^y\ ( 1 ) 

It can be seen in Eq. Q that the potential is in fact composed 
of two harmonic oscillators that has been coupled by the 
perturbation terms x^y - l/3y^. 

The Henon-Heiles potential along with the dihedral D4 
potential Q and the Toda potential l(37]| belong to a special¬ 
ized category of potentials. These potentials were obviously 
not chosen with the aim of modeling any particular type of 
galaxy. Instead, in the pioneering spirit of Henon & Heiles, 
the aim was to select potentials which are generic in their ba¬ 
sic properties, but convenient computationally, so that large 
numbers of computations could be performed. Furthermore, 
the Henon-Heiles potential admits a discrete triangular rota¬ 
tion D3 (27113) symmetry. 

The basic equations of motion for a test particle with a 
unit mass are 

V = - -— = - V - 2xy, 
ox 

dV o o 

y = -^ = -^-y + y, ( 2 ) 

dy 

where, as usual, the dot indicates derivative with respect to 
the time. Furthermore, the variational equations governing 
the evolution of a deviation vecto][^w = (dx, 6y, Sx, Sy) are 

(dv) = Sx, 

(Sy) = Sy, 

d^V d^V 

(Sx) = -^Sx - ^-^Sy = - (1 + 2y) Sx - 2xSy, 
ox^ oxoy 

(Sf) = -^-^Sx - -^Sy = -2xSx - (1 - 2y)Sy. (3) 
oyox oy^ 

Consequently, the Hamiltonian to potential 0 reads 

H - i -r -r -r -r x^y - iy^ = h, (4) 

where x and y are the momenta per unit mass conjugate to x 
and y, respectively, while /z > 0 is the numerical value of the 
Hamiltonian, which is conserved. 

The Henon-Heiles potential has a finite energy of escape 
(hesc) which is equal to 1/6. For values of energy h < hesc^ 
the equipotential curves of the system are close thus making 
escape impossible. However, for larger energy levels (h > 
hesc), the equipotential curves open and three exit channels 
appear through which the test particles may escape to infin¬ 
ity. The equipotential curves of the Henon-Heiles potential 
for various values of the energy h are shown in Fig. [^. In 

^ If is the 2N dimensional phase space (four dimensional, in our 
case) where the orbits of a dynamical system evolve on, then a devia¬ 
tion vector w, which describes a small perturbation of a specific orbit 
X, evolves on a 2N dimensional space T^S tangent to S . 
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Fig. 1 (a-upper left): Equipotential curves of the Henon-Heiles potential for various values of the energy h. The equipotential curve corresponding 
to the energy of escape is shown with red color, while the positions of the three saddle points are marked by blue dots; (b-upper right): The 3D 
surface of the potential V{x,y); (c-lower left): The potential as a function of x when y = 0; (d-lower right): The potential as a function of y when 

X = 0. 


the same plot, the equipotential corresponding to the energy 
of escape hesc is plotted with red color. 


This potential has also some additional very interesting 
features. A three-dimensional surface of the potential V(x, y) 
is presented in Fig.[^. For a constant value of y, V(x, y) has 
the form of a parabola, V{x,y) = (y + 0.5)x^ -f ^(y) (the 
parabola when y = 0 is plotted in Fig.[T]:, while the potential 
as a function of y when x = 0 is shown in Fig.[^). Further¬ 
more, the potential has a stable equilibrium point at (x, y) = 
(0,0) and three saddle points: (x,y) = {(0,1), (- V3/2,1/2), 
(V3/2,1/2)}. The saddle points constitute the three corners 
of the equipotential curve V(x,y) = 1/6, that can be seen in 
Fig. [^. This triangular area that is bounded by the equipo¬ 
tential curve with energy h = 1/6 is the area of interest in 
our research. 


3 Description of the computational methods 

In order to explore the orbital structure of the Henon-Heiles 
system, we need to define samples of orbits whose proper¬ 
ties (regularity or chaos) will be identified. For this purpose, 
we define for each energy level (all tested energy levels are 
below the escape energy), dense grids of initial conditions 
regularly distributed in the area allowed by the value of the 
energy. Our investigation takes place both in the physical 
(x, y) and the phase (y, y) space for a better understanding 
of the orbital structure. In both cases, the step separation of 
the initial conditions along the x, y and y and y axes, respec¬ 
tively (in other words the density of the grids) was controlled 
in such a way that always there are at least 50000 orbits 
in every grid. For each initial condition, we integrated the 
equations of motion 0 as well as the variational equations 
^ using a double precision Bulirsch-Stoer FORTRAN al¬ 
gorithm (e.g., 1^ ) with a small time step of order of 10“^, 
which is sufficient enough for the desired accuracy of our 
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computations (i.e., our results practically do not change by 
halving the time step). In all cases, the energy integral (Eq. 
(|^) was conserved better than one part in 10“^^ although 
for most orbits it was better than one part in 10“^^. 

Knowing whether an orbit is regular or chaotic is an is¬ 
sue of paramount importance. Over the years, several chaos 
indicators have been developed in order to determine the 
nature of orbits. In our investigation, we chose to use the 
Smaller ALingment Index (SALI) method 13^ . The SALI 
has been proved a very fast, reliable and effective tool, which 
is defined as 

SALI(t) = min{d_,d+}, (5) 

where = ||wi(0 - W2(0II and d+ = ||wi(0 -r W2(0II are the 
alignments indices, while Wi(0 and W2(0, are two deviation 
vectors which initially point in two random directions. For 
distinguishing between ordered and chaotic motion, all we 
have to do is to compute the SALI for a relatively short time 
interval of numerical integration tynax- In particular, we track 
simultaneously the time-evolution of the main orbit itself as 
well as the two deviation vectors wi(0 and W2(0 in order to 
compute the SALI. The variational equations ([^, as usual, 
are used for the evolution and computation of the deviation 
vectors. We should point out that the main feature of SALI 
which essentially distinguishes it from the computation of 
the maximal Lyapunov characteristic exponent (mLCE) is 
the introduction of one additional deviation vector with re¬ 
spect to a reference orbit. Indeed, by considering the relation 
between two deviation vectors (instead of one deviation vec¬ 
tor and the reference orbit in the case of the mLCE), one is 
able to circumvent the difficulty of the slow convergence of 
mLCE to zero (or non-zero) values as time tends to infinity. 
In other word, SALI is computationally faster than mLCE, 
which requires longer integration times in order to it reveal 
the true chaotic nature of an orbit (see e.g., ED). 

The time-evolution of SALI strongly depends on the na¬ 
ture of the computed orbit since when the orbit is regular 
the SALI exhibits small fluctuations around non zero val¬ 
ues, while on the other hand, in the case of chaotic orbits the 
SALI after a small transient period it tends exponentially to 
zero approaching the limit of the accuracy of the computer 
(10“^^). Therefore, the particular time-evolution of the SALI 
allow us to distinguish fast and safely between regular and 
chaotic motion. The time-evolution of a regular (R) and a 
chaotic (C) orbit for a time period of 10"^ time units is pre¬ 
sented in Fig. We observe, that both regular and chaotic 
orbits exhibit the expected behavior. Nevertheless, we have 
to define a specific numerical threshold value for determin¬ 
ing the transition from regularity to chaos. After conducting 
extensive numerical experiments, integrating many sets of 
orbits, we conclude that a safe threshold value for the SALI 
taking into account the total integration time of 10"^ time 
units is the value 10“^. The horizontal, blue, dashed line in 



Fig. 2 Time-evolution of the SALI of a regular orbit (green color - 
R), a sticky orbit (orange color - S) and a chaotic orbit (red color - C) 
in the Henon-Heiles system for a time period of 10"^ time units. The 
horizontal, blue, dashed line corresponds to the threshold value 10“^ 
which separates regular from chaotic motion. The chaotic orbit needs 
only about 150 time units in order to cross the threshold, while on the 
other hand, the sticky orbit requires a considerable longer integration 
time of about 6100 time units so as to reveal its true chaotic nature. 



-0.4 -0.2 0.0 0.2 0.4 0.6 


y 

Fig. 3 Regions of different values of the SALI on the (y, y) phase plane 
when the orbital energy is /z = 0.12. Light reddish colors correspond to 
ordered motion, dark blue/purplr colors indicate chaotic motion, while 
all intermediate colors suggest sticky orbits. The three black dots point 
the initial conditions of the orbits of Eig.[^ 

Fig. corresponds to that threshold value which separates 
regular from chaotic motion. In order to decide whether an 
orbit is regular or chaotic, one may use the usual method 
according to which we check after a certain and predefined 
time interval of numerical integration, if the value of SALI 
has become less than the established threshold value. There- 
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fore, if SALI < 10“^ the orbit is chaotic, while if SALI 
> 10"^ the orbit is regular. In Fig. we present a dense grid 
of initial conditions (yo^yo)^ where the values of the SALI 
are plotted using different colors. We clearly observe sev¬ 
eral regions of regularity indicated by light reddish colors 
as well as a unified chaotic domain (blue/purple dots). All 
intermediate colors correspond to sticky orbits. The initial 
conditions of the three orbits (regular, sticky and chaotic) 
whose SALI time-evolution was given in Fig. are pin¬ 
pointed with black dots in Fig.[^ Therefore, the distinction 
between regular and chaotic motion is clear and beyond any 
doubt when using the SALI method. 

In our study, each orbit was integrated numerically for 
a time interval of 10 "^ time units (10 billion yr), which cor¬ 
responds to a time span of the order of hundreds of orbital 
periods. The particular choice of the total integration time 
is an element of great importance, especially in the case of 
the so called “sticky orbits” (i.e., chaotic orbits that behave 
as regular ones during long periods of time). A sticky orbit 
could be easily misclassified as regular by any chaos indi- 
catoij^ if the total integration interval is too small, so that 
the orbit do not have enough time in order to reveal its true 
chaotic character. Thus, all the sets of orbits of a given grid 
were integrated, as we already said, for lO'^ time units, thus 
avoiding sticky orbits with a stickiness at least of the order 
of a Hubble time. All the sticky orbits which do not show 
any signs of chaoticity for 10^ time units are counted as reg¬ 
ular ones, since that vast sticky periods are completely out of 
scope of our research. A characteristic example of a sticky 
orbit (S) in our galactic model can be seen in Fig. where 
we observe that the chaotic character of the particular sticky 
orbit is revealed only after a considerable long integration 
time of about 6100 time units. 

A first step towards the understanding of the overall or¬ 
bital behavior of our galactic system is knowing the regular 
or chaotic nature of orbits. In addition, of particular inter¬ 
est is the distribution of regular orbits into different families. 
Therefore, once the orbits have been characterized as regular 
or chaotic, we then further classified the regular orbits into 
different families, by using the frequency analysis of ca. 
Initially, (TOl [TTl proposed a technique, dubbed spectral dy¬ 
namics, for this particular purpose. Later on, this method has 
been extended and improved by Ea and oa. In a recent 
work, l42l the algorithm was refined even further so it can 
be used to classify orbits in the meridional plane. In general 
terms, this method calculates the Fourier transform of the 
coordinates of an orbit, identifies its peaks, extracts the cor¬ 
responding frequencies and search for the fundamental fre- 


^ Generally, dynamical methods are broadly split into two types: (i) 
those based on the evolution of sets of deviation vectors in order to 
characterize an orbit and (ii) those based on the frequencies of the or¬ 
bits which extract information about the nature of motion only through 
the basic orbital elements without the use of deviation vectors. 


Table 1 Type, initial conditions and the value of the energy of the or¬ 
bits shown in Fig. [4la-h). In all cases, xq = 0, xq is found from the 
energy integral (EqTQ), while Tper is the period of the resonant parent 
periodic orbits. 


Figure Type 

h 

yo 

yb 

Tper 

4 

i 

box 

0.02 

-0.08200000 0.00000000 

- 

4 


1:1 loop 

0.12 

0.29754452 

0.00000000 

6.08239553 

4 


1:1 linear 

0.12 

0.00000000 

0.24494897 

7.43702045 

4 

1 

2:2 

0.16 

0.50003040 

0.00000000 

12.47212865 

4 


3:3 

0.14 

0.49779098 

0.00000000 

18.55933975 

4 


4:4 

0.12 

0.10124417 

0.28443422 

29.16285964 

4 


5:5 

0.12 

0.54846395 

0.00000000 

32.13261997 

4 


chaotic 

0.16 

0.82000000 

0.00000000 

- 


quencies and their possible resonances. Thus, we can easily 
identify the various families of regular orbits and also recog¬ 
nize the secondary resonances that bifurcate from them. The 
very same algorithm was used in several previous related 
papers (e.g., ll4^rT5ll40ll4T]l L 

4 Numerical results & Orbit classification 

The main objective of this numerical investigation is to dis¬ 
tinguish between ordered and chaotic orbits for values of 
energy always lower than the escape energy. In order to ac¬ 
complish a spherical view of the orbital properties of the 
Hamiltonian system, we expand our investigation in both 
the physical and the phase space. In both cases, the grids of 
initial conditions of orbits whose properties will be exam¬ 
ined are defined as follows: for the physical (v,y) space we 
consider orbits with initial conditions (vo,yo) with yh = 0. 
while the initial value of xq is obtained from the energy in¬ 
tegral 0 as Xq = x(xo,yo,h) > 0. Similarly, for the phase 
(y,y) space we consider orbits with initial conditions (yo,yo) 
with Xq = 0, while again the initial value of xq is obtained 
from the energy integral 0. As we explained earlier, in 
our investigation we shall deal only with bounded motion 
(h e (0, hesc)) of tost particles for values of energy in the set 
h = {0.02,0.04,0.06,0.08,0.10,0.12,0.14,0.16}. The case 
of unbounded motion (h > hesc) iii the Henon-Heiles system 
has been investigated extensively in many previous papers 
(e.g., 

At this point, we would like to clarify some issues re¬ 
garding the nomenclature of the orbits in our model. In the 
present paper, we decided to follow the classification method 
according to which, the orbits are separated into three main 
categories: (i) box orbits, (ii) n : m resonant orbits, and (iii) 
chaotic orbits. According to our notation, all resonant orbits 
have the following recognizable n : m oscillatory pattern: a 
resonant orbit completes m oscillations perpendicular to the 
X axis in the time that it takes the orbit to perform n circuits 
along the y axis. Furthermore, slh : m resonant orbit would 
be represented by n distinct islands of invariant curves in the 
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Fig. 4 Collection of the basic types of orbits encountered in the Henon-Heiles system: (a) box orbit; (b) 1:1 resonant loop orbit; (c) 1:1 resonant 
linear orbit; (d) 2:2 resonant orbit; (e) 3:3 resonant orbit (f) 4:4 resonant orbit; (g) 5:5 resonant orbit; (h) chaotic orbit. 




h h 

Fig. 5 (a-left): The (yo, h) characteristic curves of the orbital families and (b-right): Evolution of the stability index S.L of the families of periodic 
orbits shown in Eigs.|^. The black dashed horizontal lines at -2 and +2 delimit the range of S.L for which the periodic orbits are stable. The color 
code is the same in all panels. 


(x, x) phase plane and m distinct islands of invariant curves 
in the (y,y) surface of section (e.g., ifT^ ). The n : m nota¬ 
tion we use for the regular orbits is according to (Tbl l42ll . 
where the ratio of those integers corresponds to the ratio of 
the main frequencies of the orbit, where main frequency is 
the frequency of greatest amplitude in each coordinate. Main 
amplitudes, when having a rational ratio, define the reso¬ 
nances of an orbit. In our research, we searched for resonant 
orbits n : mup to n,m < 10 therefore, for all higher resonant 
orbits (if any) the numerical code assigns “box” classifica¬ 
tion (this is a usual technique in orbit classification), which 


is indeed correct ioxn ^ m (high resonant box orbits (e.g., 

d). 

Our numerical calculations suggest that in the Henon- 
Heiles system apart from the usual box and chaotic orbits 
there are also many types of resonant periodic orbits. The 
majority of resonant orbits belong to the main 1:1 resonant 
family however, several secondary bifurcated families (i.e., 
the 2:2, 3:3, 4:4 and 5:5) are also present. Moreover, it was 
observed that the contribution of other resonant orbits (i.e., 
2:1, 3:2, or higher) to the overall orbital structure of the sys¬ 
tem is extremely limited. In Fig.j^a-h) we present an exam- 
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Fig. 6 The orbital structure of the physical (x, y) plane for several values of the energy h, distinguishing between chaotic and different types of 
regular orbits. 



Fig. 7 The orbital structure of the phase (y,y) plane for several values of the energy h, distinguishing between chaotic and different types of regular 
orbits. 


pie of each of the basic types of the orbits. The box orbit and 
the chaotic orbit shown in Figs. Eh an d|4}i, respectively were 
computed until t = 400 time units, while all the parenj^ pe¬ 
riodic orbits were computed until one period has completed. 
The black thick curve circumscribing each orbit in the phys¬ 
ical (v,y) plane is the Zero Velocity Curve (ZVC) which is 


^ For every orbital family there is a parent (or mother) periodic or¬ 
bit, that is, an orbit that describes a closed figure. Perturbing the initial 
conditions which define the exact position of a periodic orbit we gen¬ 
erate quasi-periodic orbits that belong to the same orbital family and 
librate around their closed parent periodic orbit. 


defined as V{x,y) = h. In Table we provide the type, the 
initial conditions and the value of the energy for all the de¬ 
picted orbits. In the resonant cases, the initial conditions and 
the period Tper correspond to the parent periodic orbits. 

A very informative diagram the so-called “characteris¬ 
tic” orbital diagram ^T2\ is presented in Fig. [^. This dia¬ 
gram shows the evolution of the y coordinate of the initial 
conditions of the parent periodic orbits of each orbital fam¬ 
ily as a function of their orbits energy h. Here we should 
emphasize, that for orbits starting perpendicular to the y- 
axis, we need only the initial condition of yo in order to lo- 
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h h 

Fig. 8 Evolution of the percentages of the different types of orbits when varying the energy h (a-left): on the physical (x,y) plane and (b-right): on 
the phase (y, y) plane. 


cate them on this characteristic diagram. On the other hand, 
for orbits not starting perpendicular to the y-axis (i.e., the 
1:1 linear and the 4:4 resonant orbits) initial conditions as 
position-velocity pairs (y, y) are required and therefore, the 
characteristic diagram is now three-dimensional providing 
the full information regarding the interrelations of the initial 
conditions in a tree of families of periodic orbits. Further¬ 
more, the diagram shown in Fig. is called the “stability 
diagram” [21 and it illustrates the stability of all the 
families of periodic orbits in our dynamical system when 
the numerical value of the energy varies. A periodic orbit is 
stable if only the stability index (S.I.) Il26ll?9l is between 2 
and -f2. This diagram help us monitor the evolution of S.I. 
of the resonant periodic orbits as well as the transitions from 
stability to instability and vice versa. 

In the following Figs. [^and[7] we present six grids of 
initial conditions (vo,yo) and (yo,yo) of orbits that we have 
classified for different values of the total orbital energy h in 
the physical and in the phase space, respectively. The corre¬ 
sponding black thick curve circumscribing each grid in the 
phase (y, y) plane is the limiting curve which can be obtained 
from 

fiy,y) = \f + V{x = 0,y) = h. (6) 

We should point out, that the grids in both the physical and 
the phase space have the form of a nut, being shaper in its 
comers when the value of energy is close to escape energy, 
while it gets smoother and more spherical as h goes to zero. 
Looking at Figs. and |7] it is evident that at very low en¬ 
ergy levels the motion is entirely ordered. In particular, the 
vast majority of the regular orbits are either box or 1:1 reso¬ 
nant orbits which form three distinct stability islands. We 
also see that with increasing energy, the total area in the 


grids occupied by box orbits is constantly reduced and for 
h > 0.06 initial conditions corresponding to box orbits are 
mainly confined to the outer parts of the grids. Moreover, 
it is seen that when h = 0.08 a thin chaotic layer emerges 
inside the vast box region. For higher energy levels, the or¬ 
bital structure of the computed grids alters significantly and 
the most noticeable changes are the following: (i) the initial 
conditions of box orbits practically vanish; (ii) the area oc¬ 
cupied by 1:1 resonant orbits is reduced; (iii) the amount of 
chaotic orbits is highly increased and (iv) several bifurcated 
families of the main 1:1 family (i.e., the 2:2, 3:3, 4:4 and 5:5 
resonant families) along with other higher resonant orbital 
families appear. Furthermore, we observe that for h > 0.12, 
that is when the value of the energy approaches very close to 
the escape energy, chaotic orbits is the most populated type 
of orbits however, small stability islands are still present and 
embedded inside the vast unified chaotic sea. 

The evolution of the percentages of the chaotic orbits 
and that of the main families of regular orbits on the physical 
as well as on the phase space when the total orbital energy 
h varies is shown in Fig.[^a-b). We observe, that the pattern 
of the curves is more or less the same in both the physical 
and the phase space. It is seen that for low values of the 
energy there is no evidence of chaotic motion. In particular, 
only two type of regular orbits are possible: box orbits which 
cover about 70% of the total initial conditions and 1:1 res¬ 
onant orbits which correspond to the remaining 30% of the 
computed orbits. With increasing energy however, the per¬ 
centage of the box drops rapidly, while at the same time the 
rate of the 1:1 resonant orbits grows significantly and very 
quickly (h > 0.03) the 1:1 resonant orbits is the most popu¬ 
lated family. In fact, in the interval 0.06 < h < 0.10, the 1:1 
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Fig. 9 The orbital structure of the (a-left): (x,/z)-plane and (b-right): (y,/z)-plane. These diagrams give a detailed analysis of the evolution of the 
regular and chaotic orbits of the dynamical system when the value of the energy h changes. 
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Fig. 10 The orbital structure of the (a-left): (x, /z)-plane and (b-right): (y, /z)-plane. These diagrams give a detailed analysis of the evolution of the 
chaotic as well as of the different types of regular orbits of the dynamical system when the value of the energy h changes. 


resonant orbits occupy more than 90% of the grids, while for 
larger values of the energy this tendency is reversed, their 
rate reduces rapidly and for /z = 0.16 1:1 resonant orbits 
occupy only about one tenth of the grids. The rate of box 
orbits on the other hand, is constantly being reduced and 
tends asymptotically to zero as we proceed to higher energy 
levels. Furthermore, we may say that chaos appears only at 
relatively high energies and when h > 0.08 the percentage 
of chaotic orbits exhibits a drastic increase and eventually 
chaotic orbits take over the grids when /z > 0. 12|^ The par¬ 
ticular evolution of the percentage of chaotic orbits was an- 

^ It is known that in dynamical systems with a finite energy of es¬ 
cape where escape is possible (open Hamiltonian systems), chaos ap¬ 
pears and dominates at values of energy very close to the critical escape 


ticipated because as the value of the energy increases we 
approach closer and closer to the critical value of the escape 
energy (hesc = 1/6). Indeed, at the higher energy value stud¬ 
ied, that is when h = 0.16, about 90% of the total initial 
conditions correspond to chaotic orbits. Moreover, we may 
argue that varying the value of the energy mainly shuffles the 
orbital content of all the other resonant orbits, whose per¬ 
centages present fluctuations at low values (always less than 
10%). However, we would like to note that when h = 0.12, 
we observe a sudden spike at the percentages of 4:4 and 5:5 
resonant orbits which belong to bifurcated families of the 
main 1:1 resonant family. In addition, the rate of the 2:2 res- 

energy which determines the transition from bounded to unbounded 
motion. 
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onant orbits is slightly elevated at /z = 0.16 thus surpassing 
that of the main 1:1 family. Thus, taking into account ah 
the above-mentioned analysis, we may conclude that in the 
Henon-Heiles system the total orbital energy h influences 
mostly box, 1:1, and chaotic orbits. In fact, what really hap¬ 
pens is that a large portion of 1:1 and box orbits turns into 
chaotic orbits as the value of the energy increases tending to 
the escape energy. 


The grids in the physical (x,y) as well as in the phase 
(y, y) plane provide information on the phase space mixing 
for only a flxed value of energy. Henon however, back in 
the 60s, also considered two different types of planes which 
provide information about regions of regularity and chaos, 
using the value of the energy as an ordinate thus obtaining 
a continuous spectrum of values of h. The first type which 
is the (v,/z) plane uses the section y = i: = 0, y > 0, i.e., 
the test particle starts on the v-axis, parallel to the y-axis and 
in the positive y-direction. Similarly, the second type, the 
(y, h) plane uses the section x = y = 0, i: > 0, i.e., the test 
particle starts on the y-axis, parallel to the x-axis and in the 
positive x-direction. In Fig. [^a-b) we show the structure of 
these planes when h e (0, hesc)^ where each initial condition 
is plotted in color code according to the value of SALI, thus 
distinguishing between ordered and chaotic motion. We ob¬ 
serve that the vast majority of the computed orbits are reg¬ 
ular, while on the other hand, initial conditions of chaotic 
orbits are mainly confined to right outer parts of the grids. 
In order to be able to monitor with sufhcient accuracy and 
details the transition from regularity to chaos, we defined a 
dense grid of 10^ initial conditions in both cases. Here we 
should point out, that our results shown in Fig. [^a-b) coin¬ 
cide with that of Figs. 4 and 6, respectively presented in lO, 
where the OFLI2 indicator was deployed in order to disso¬ 
ciate regular and chaotic orbits. 


Finally, Fig. [^a-b) presents another point of view of 
the planes discussed earlier in Fig. [^a-b). In this case, the 
initial conditions of the orbits are colored according to the 
regular or chaotic nature nature of the orbits giving empha¬ 
sis to the different types of regular orbits. Once more, we 
observe that at low energy levels the vast majority of orbits 
are regular and speciflcahy box orbits. As the value of the 
energy increases however, box orbits are depleted and 1:1 
resonant orbits is the most populated family of regular or¬ 
bits. It is also seen, that ah the bifurcated families (i.e., the 
2:2, 3:3, 4:4 and 5:5 families), as well as other higher res¬ 
onant families of orbits (i.e., the 3:2 family) emerge at rel¬ 
atively high values of the energy (roughly about h > 0.08). 
Around the same energy level we also encounter the first 
initial conditions corresponding to chaotic motion. The hor¬ 
izontal black dashed line at Fig. [^a-b) marks the first in¬ 
dication of chaos. Our numerical calculations suggest that 
at the (x,/z) plane chaotic motion appears at /z = 0.07916, 


while at the (y, h) plane the energy level in which chaos was 
detected for the first time equals to /z = 0.07882. 


5 Discussion and conclusions 

The aim of the present work was to investigate how influen¬ 
tial is the total orbital energy on the level of chaos and on 
the distribution of the different regular families in the classi¬ 
cal Henon-Heiles system. A very interesting feature of this 
Hamiltonian system is that it has a finite energy of escape. 
In particular, for energies smaller than the escape value, the 
equipotential surfaces are close and therefore escape is im¬ 
possible. For energy levels larger than the escape energy on 
the other hand, the equipotential surfaces open and three 
channels of escape appear through which the particles can 
escape to infinity. However, it should be emphasized, that if 
a test particle has energy larger than the escape value, this 
does not necessarily mean that the particle will certainly es¬ 
cape from the system. In our research, we decided to con¬ 
sider only bounded motion and explore in details the regular 
or chaotic nature of orbits. 

We numerically investigated the structure of both the 
physical (x, y) and the phase (y, y) space for a better under¬ 
standing of the orbital properties of the system. For deter¬ 
mining the regular or chaotic character of orbits in this sys¬ 
tem, we defined dense grids of initial conditions regularly 
distributed in the area allowed by the value of the total or¬ 
bital energy /z, in both the physical and the phase space. In 
both cases, the density of the grids was controlled in such a 
way that always there are at least 50000 orbits to be exam¬ 
ined. For the numerical integration of the total set of orbits 
in each grid, we needed roughly between 17 to 40 hours of 
CPU time on a Pentium Dual-Core 2.2 GHz PC, depending 
mainly on the amount of regular orbits in each case which, 
unlike chaotic orbits, they required further numerical calcu¬ 
lations in order to be classified into regular families. 

To show how the value of the total energy h influences 
the orbital structure of the system, we presented for each 
case, color-coded grids of initial conditions, which allow us 
to visualize what types of orbits occupy specific areas in the 
physical/phase space. Each orbit was integrated numerically 
for a time period of 10"^ time units which corresponds to a 
time span of the order of hundreds of orbital periods. The 
particular choice of the total integration time was made in 
order to eliminate sticky orbits (classifying them correctly 
as chaotic orbits) with relatively long sticky periods. Then, 
we made a step further, in an attempt to distribute ah regu¬ 
lar orbits into different families. Therefore, once an orbit has 
been characterized as regular applying the SALI method, we 
then further classified it using a frequency analysis method. 
This method calculates the Fourier transform of the coordi¬ 
nates and velocities of an orbit, identifies its peaks, extracts 
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the corresponding frequencies and then searches for the fun¬ 
damental frequencies and their possible resonances. 

Our thorough and detailed numerical analysis suggests 
that the level of chaos as well as the distribution in regular 
families are indeed very dependent on the total orbital en¬ 
ergy. The main numerical results of our investigation can be 
summarized as follows: 

1. Several types of regular orbits were found to exist in 
Henon-Heiles potential, while there are also extended 
chaotic domains separating the areas of regularity. In 
particular, a large variety of resonant orbits are present, 
thus making the orbital structure of the system both rich 
and interesting. 

2. Our numerical computations revealed that most of the 
resonant families (i.e., the 2:2, 3:3, 4:4, 5:5 families) are 
in fact bifurcated families of the main 1:1 resonant fam¬ 
ily. Other types of resonant orbits on the other hand, (i.e., 
the 2:1, 3:2, etc) have extremely weak contribution to the 
overall orbital structure of the system. 

3. It was found that in both the physical (x, y) and the phase 
(y,y) space the total orbital energy influences mainly 
box, 1:1 resonant, and chaotic orbits. Moreover, for h > 
0.06 the majority of test particles move either in 1:1 res¬ 
onant or chaotic orbits. 

4. For low energy levels (h < 0.08) the motion was found to 
be entirely regular, while for larger values of the energy 
the amount of chaotic orbits increases rapidly and for 
h> 0.14 they dominate occupying more than half of the 
entire phase space. 

5. The (x, h) and (y, h) diagrams gave us a more complete 
view of the evolution of the orbital families for a con¬ 
tinuous spectrum of values of h. It was observed, that 
for low energy levels corresponding to local motion the 
orbital content is rather poor (only box and 1:1 reso¬ 
nant orbits are present), while most of the resonant or¬ 
bits emerge at relatively high values of energy (/t > 0.08) 
suitable for global motion. 

We hope that the present analysis and the corresponding 
numerical results of the Henon-Heiles system to be useful 
in this active held of nonlinear Hamiltonian systems. Taking 
into account that our outcomes are encouraging, it is also in 
our future plans to expand our exploration in order to reveal 
the complete network of periodic orbits, thus shedding some 
light to the evolution of the periodic orbits as well as their 
stability when the orbital energy varies, considering also un¬ 
bounded motion of particles. 
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